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Abstract. The results of analytical approximations and extensive cal- 
culations based on a path integral Monte Carlo (PIMC) scheme are pre- 
sented. A new (direct) PIMC method allows for a correct determination 
of thermodynamic properties such as energy and equation of state of 
dense degenerate Coulomb systems. In this paper, we present results 
for dense partially ionized hydrogen at intermediate and high temper- 
ature. We give a quantitative comparison with the available results of 
alternative (restricted) PIMC simulations and with analytical expres- 
sions based on iterpolation formulas meeting the exact limits at low and 
high densities. Good agreement between the two simulations is found up 
to densities of the order of 10 24 cm -3 . The agreement with the analytical 
results is satisfactory up to densities in the range 10 22 . 



1 Introduction 

Correlated Fermi systems are of increasing interest in many fields, including 
plasmas, astrophysics, solids and nuclear matter, (see Kraeft et al. 1986) for an 
overview. Among the topics of current interest are Fermi liquids, metallic hy- 
drogen (see DaSilva et al. 1997), plasma phase transition (see Schlanges et al. 
1995), bound states etc. In such many particle quantum systems, the Coulomb 
interaction is essential. There has been significant progress in recent years to 
study these systems theoretically, and especially numerically, (see e.g. Bonitz 
(Ed.) 2000, Zamalin et al. 1977, Filinov, A. V. et al. 2000). A theoretical frame- 
work which is particularly well suited to describe thermodynamic properties 
in the region of strong coupling and degeneracy is the path integral quantum 
Monte Carlo (PIMC) method. There has been remarkable recent progress in 
applying these techniques to Fermi systems. However, these simulations are es- 
sentially hampered by the fermion sign problem. To overcome this difficulty, 
several strategies have been developed to simulate macroscopic Coulomb sys- 
tems (see Militzer and Pollock 2000, Militzer and Ceperley 2000, and Militzer 
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2000) : the first is the restricted PIMC concept where additional assumptions on 
the density operator p are introduced which reduce the sum over permutations 
to even (positive) contributions only. This requires the knowledge of the nodes 
of the density matrix which is available only in a few special cases. However, 
for interacting macroscopic systems, these nodes arc known only approximately, 
(see, e.g., Militzer and Pollock 2000 and Militzcr and Ccpcrlcy 2000), and the 
accuracy of the results is difficult to assess from within this scheme. 

Recently, we have published a new path integral representation for the N- 
particle density operator (see Filinov, V. S., et al. 2000, Bonitz (Ed.) 2000, 
Filinov, V. S., et al. 2001), which allows for direct Fermionic path integral Monte 
Carlo simulations of dense plasmas in a broad range of densities and tempera- 
tures. Using this concept we computed the pressure (equation of state, EOS), 
the energy, and the pair distribution functions of a dense partially ionized and 
dissociated electron-proten plasma (see Filinov, V. S., et al. 2000). In this re- 
gion no reliable data are available from other theories such as density functional 
theory or quantum statistics (see, e.g., Kraeft et al. 1986), which would allow 
for an unambiguous test. Therefore, it is of high interest to perform quantitative 
comparisons of analytical results and independent numerical simulations, such 
as restricted and direct fermionic PIMC, which is the aim of this paper. 



2 Path integral representation of thermodynamic 
quantities 

We now briefly outline the idea of our direct PIMC scheme. All thermodynamic 
properties of a two-component plasma are defined by the partition function Z 
which, for the case of N e electrons and N p protons, is given by 



N e \N p \ 



with Q(N e ,N p ,(3) =J2 J dqdrp(q,r,a;(3) 7 



(1) 



where [3 — 1/ksT. The exact density matrix is, for a quantum system, in general, 
not known but can be constructed using a path integral representation (see 
Feynman and Hibbs 1965), 

J dRM ]T p(i? (0) , a; 0) = J dRM . . . dR™ p (1) • P {2) . . . P {n) 
v v 

x££(±i)^s(*,A/)fy»+ i > ) (2) 

a P 

where p« = p (i^" 1 ), A0) = {R^-^\e- A ^\R^}, whereas Aj3 = f3/(n + 
1). H is the Hamilton operator, H = K + U Cl containing kinetic and potential en- 
ergy contributions with U c = lJ^ + Uf + JJ^ p being the sum of the Coulomb poten- 
tials between protons (p), electrons (e) and electrons and protons (ep). Further, 
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R (i) = (g« r «) = for i= l,...n+l, = for) = (R P 0) , R { ° ] ). 



Also, R( n+1 ) = R(°> and a' = a, i.e., the particles are represented by closed 
Fermionic loops with the coordinates (beads) [R] = [RW;RW; . . . ; R^; R ( - n+1 '>}, 
where r and g denote the electron and proton coordinates, respectively. The spin 
gives rise to the spin part of the density matrix S, whereas exchange effects are 
accounted for by the permutation operator P, which acts on the electron co- 
ordinates and spin, and the sum over the permutations with parity up. In the 
fermionic case (minus sign), the sum contains N e \/2 positive and negative terms 
leading to the notorious sign problem. Due to the large mass difference of elec- 
trons and protons, the exchange of the latter is not included. 

Recently, we have derived a new representation for the high-temperature 
density matrices pw in eq.(2) (see Filinov, V. S., et al. 2000) which is well 
suited for direct PIMC simulations. A crucial point is that the electron-proton 
interaction can be described by an (effective) quantum pair potential <P ep {Kelbg- 
potential, see Kelbg 1964). For details of the derivation see Filinov, V. S., et al. 
2000. Here, we present only the final result for the energy and for the EOS. 
Consider first the energy: 
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and <Pf = A/3d[f3'$ ep (\x l pt \,/3')]/d/3'\ >=Af3 contains the electron-proton Kelbg 
potential <P ep . Here, (. . . | . . .) denotes the scalar product, and q pt , r pt and x p t are 
differences of two coordinate vectors: q p t = q p — qt, r pt = r p — r t , x pt = r p — q t , 
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yl with y% = A\ e Y, n k =id k) and 



AX 2 a = 2-K% z Ap/m a . We introduced dimensionless distances between neighbor- 
ing vertices on the loop, ^\ . . .S} n \ thus, explicitly, [r] = [r; y^ 1 ; ; •••]• 
The density matrix p s is given by 



n N B 



P.(?,[r],/3)=CS r .e-^W.«Jin 



^det|^ a Yl 



(4) 



i=i p=i 



where U(q, [r],/3) = U p (q) + {U e ([r}, Ap) + U ep (q, [r], Ap)}/(n + 1) and 0! 
exp[— 7r|^p | 2 ]. We underline that the density matrix (|j) does not contain an 
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explicit sum over the permutations and thus no sum of terms with alternating 
sign. Instead, the whole exchange problem is contained in a single exchange 
matrix given by 
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As a result of the spin summation, the matrix carries a subscript s denoting 
the number of electrons having the same spin projection. For more detail, (see 
Filinov, V. S., et al. 2000, Bonitz (Ed.) 2000). In a similar way, we obtain the 
result for the equation of state, 
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3 Analytical approximations for the thermodynamic 
functions of dense plasmas 

To describe dense plasmas, it is necessary to have thermodynamic functions valid 
at arbitrary degeneracy. Here, we restrict ourselves to the Hartree-Fock (HF) 
and the Montroll-Ward (MW) contributions. This approximation is appropriate 
at temperatures high enough such that the Coulomb interaction is weak and the 
possibility of the formation of bound states is excluded. HF and MW contribu- 
tions have been computed numerically (see Kraeft ct al. 1986). The analytical 
evaluation of the MW contribution is possible in limiting situations only, namely 
in the low and very high density cases. In the intermediate region Pade formulae 
can be used to interpolate between the limiting cases. In between, the formulae 
are fitted to numerical data; (see Ebeling et al. 1981 , Haronska et al. 1987 and 
Ebeling and Richert 1985). We give the excess free energy and the interaction 
part of the chemical potential of the electron gas, 

£P f D -l(7rP)-^n + 8n 2 f GB 



1 + 8ln 



64 V2 



+ 8n 2 



V 



and 



1 + 8ln 



1(^)1/4^1/2 



16V2 



8n 2 



(8) 



In (0J|) Heaviside units h = ^- = 2m e = 1 and the dimensionless density 
fi = were used. 
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Fig. 1. Energy and pressure isotherms for 50, OOOi-T (solid line). Solid line with 
circles- Hartree Fock (HF) and Montroll-Ward (MW) approximation. Reference 
data: triangle - RPIMC (see Militzer et al.), square - DPIMC. 



In formulae (0,^), the correct low density behaviour (Debye limiting law) is 
guaranteed by choosing fp — — (2/3)(7r/3) _1 / 4 n 1 / 2 and fiu = —{■K(3)~ 1 / i n 1 / 2 . 
The correct high degeneracy limit is recovered by using the (slightly modified) 
Gell-Mann Brueckner approximations (including Hartree-Fock) for the free en- 
ergy and for the chemical potential 
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The free energy is now equal to the internal energy atT = and reads, according 
to Carr and Maradudin, 

U 2 21 916 

— = — -= : + 0.06221nr s - 0.096 + 0.018r s lnr s + ■ ■ • . 

N rj r s 

The Brueckner parameter r s is given by rf = 3/(47rn). While the Hartree Fock 
term, i.e., the l/r s term in (|9|,[l0|), was retained unaffected, the additional terms 
in these equations and in formulae (0,^) were modified, or fitted, respectively, 
such that (|7],||) meet the numerical data in between, where the analytical limiting 
formulae are not applicable. 
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Fig. 2. Energy and pressure isotherms for 100, 000K (solid line with circles). 
Solid line with small squares - Hartree Fock (HF) and Montroll-Ward (MW) 
approximation. Reference data: triangle - RPIMC (see Militzer et al. 2000), large 
square - DPIMC. 
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Consider now the proton contributions. In the low density regime, the proton 
formulae are practically the same as for the electrons, whereas in the high density 
limit we adjust the formulae to (classical) Monte-Carlo (MC) data. We have for 
the free energy density and for the chemical potential 

U _ (-f p nt /k B Tn p ) D [l - an 2 p /3 (f^/k B Tn p )Mc} 



fi p _ (-»f/k B T) D [l - 2an 2 p / \^/k B Tn p ) MC ] 

knT ~ i o ~l/2 r -l/2 //l4 nt \ . -1/6/ACS 1 

We used the abbreviation a which depends only on temperature 
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a=V^W\l\l + \ '^-e X] >(— ^ _)1 - 0.29931}. (13) 

For the protons, we introduce the dimensionless density to be used in the plasma 
parameter r, namely n p = n^pms n p , r — ( ^Trh p ) . The Debye approxima- 
tions (i.e., the low density case) for free energy density and chemical potential 
read {- f p nt /k B Tn p ) D = 2.1605rip /2 and {-^ p nt /k B T) D = §2.1605ftp /2 . In the 
high density region, we use a fit to (classical OCP) Monte-Carlo data. 
For the free energy we write 
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and for the chemical potential 

(n p nt /k B T) MC = -1.1928P + 3.5382/ 11 / 4 - 0.5012 lnP - 2.9761 
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The correlation part of the pressure for an iJ-plasma is then given by 

p corr = p corr + ^corr = ^ + ^ _ /g _ /p . (Jg) 

The contributions are determined by (@,|) and QJH)- The ideal pressure is 
given by Fermi integrals I v {a) 

p ld = k B T^2 h/2{a a ) , a a = Ha/{k B T) . 

a a 

The internal energy may be constructed from the excess free energy / = y given 
above in addition to the ideal part according to U — F — ^fff|y_ const > where 
the ideal free energy is given by 
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At very high degeneracy, the free energy is equal to the internal energy. 



4 Hydrogen isotherms 

In this section we present results for the thermodynamic functions of dense 
hydrogen versus density at constant temperature. The PIMC simulations have 
been performed as explained in Filinov, V. S., et al. 2000, and Filinov V. S., et 
al. 2001. Figures 1-3 show the simulation results together with the Pade results 
for three hydrogen isotherms T = 50, 000; 100, 000; and 125, 000K. In all figures 
the agreement between numerical and analytical data is good for temperatures 
and densities, where the coupling parameter r is smaller than or equal to unity. 
Reference points related to RPIMC and DPIMC calculations in Fig. 1 correspond 
to data available for the temperature T — 62, 500K. 




Density n, cm' 

Fig. 3. Energy and pressure isotherms for 125, 000K (lower curves, solid line 
with larger or without circles). Solid line with smaller circles - Hartrec Fock 
(HF) and Montroll-Ward (MW) approximation for 125000.fr. Reference data: 
triangles - RPIMC (see Militzer et al. 2000). 
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At low densities, pressure and energy are close to those of an ideal plasma. 
Increasing the density above 10 19 cm~ 3 , Coulomb interaction becomes impor- 
tant leading to a decrease of pressure and energy. Differences between analytical 
and numerical calculations in Fig. 1 are observed for densities above 10 22 cm~ 3 
where the coupling parameter r exceeds unity. At temperatures of 100, 000 and 
125, OOOif, differences are observed for n above 5x 10 22 cto~ 3 . The degeneracy pa- 
rameter n e A 3 reaches here values of 0.4. At higher densities (around 10 24 cm~ 3 ) 
the degeneracy n e X 3 becomes larger than unity, and the interaction parts of 
pressure and energy decrease as compared to the respective ideal contributions, 
which leads to an increase of pressure and energy. At lower temperatures, this 
tendency is accompanied by the vanishing of bound states, i.e., a transition from 
a partially ionized plasma to a metal-like state. This tendency is correctly repro- 
duced by all methods, however the density values of this increase vary. In Fig. 
3 we compare our results with data from RPIMC simulations (see Militzer et al 
2000). Obviously, the agreement is very good up to densities below 10 24 cm~ 3 . 

5 Discussion 

This work is devoted to a Quantum Monte Carlo study of a correlated proton- 
clcctron system with degenerate electrons. We compared our direct PIMC sim- 
ulations with independent restricted PIMC results of Militzer and Ceperley and 
analytical formulae for isotherms corresponding to T = 50,000; 100,000; and 
125, 000if. The values of r and n e /l 3 are varying in a wide range of values. This 
region is of particular interest as here pressure and temperature ionization occur 
and, therefore, an accurate and consistent treatment of scattering and bound 
states is crucial. We found that the results agree sufficiently well for coupling 
parameters smaller or equal to unity. This is remarkable because analytical for- 
mulae, the DPIMC and RPIMC simulations are completely independent and use 
essentially different approximations. We, therefore, expect that these results for 
hydrogen are reliable which is the main result of the present paper. We hope 
that our simulation results allow us to derive and test improved analytical ap- 
proximations in the future. 

6 Acknowledgements 

We acknowledge support by the Deutsche Forschungsgemeinschaft (Sonderfor- 
schungsbereich 198: M.B., D.K., and W.D.K; Mercator-Programm: V.S.F.). Our 
thanks are due to W. Ebeling and M. Schlanges for stimulating discussions. 

References 

Kraeft, W.D., Kremp, D., Ebeling, W., Ropke, G.: Quantum Statistics of Charged 

Particle Systems, Akademie-Verlag Berlin and Plenum New York 1986 
Da Silva, LB. et al.,: Phys. Rev. Lett. 78 (1997) 783 



X 



Schlanges, M., Bonitz, M., Tschtschjan, A.: Contrib. Plasma Phys. 35 (1995)109 
Bonitz, M., (Ed.): Progress in Nonequlibrium Green's functions, World Scientific, Sin- 
gapore 2000 

Zamalin, V.M., Norman, G. E., Filinov, V. S.: The Monte Carlo Method in Statistical 

Thermodynamics, Nauka, Moscow 1977 (in Russian). 
Filinov, A. V., Lozovik, Yu. E., Bonitz, M.: phys. stat. sol. (b) 221 (20 00) 231; Fil- 



inov, A. V., Bonitz, M., Lozovik, Yu. E.: Phys. Rev. Lett. (2001), 



mat/0012265 



arXiv:cond- 



Militzer, B., Pollock, E.: Phys. Rev. E 61 (2000) 3470 
Militzer, B., Ceperley, D.: Phys. Rev. Lett. 85 (2000) 1890 
Militzer, B., PhD-Thesis, University of Illinois (2000) 

Filinov, V. S., Bonitz, M., Fortov, V. E.: JETP Letters, 72 (2000) 245; Filinov, V. S., 
Fortov, V. E., Bonitz, M., Kremp, D.: Phys. Lett. A 274 (2000) 228 

Filinov, V. S., Bonitz, M., Ebeling, W., Fortov, V. E.: Plasma Phys. Contr. Fusion 
(2001) 

Feynman, R. P., Hibbs, A. R.: Quantum mechanics and path integrals, 1965, McGraw- 
Hill, New York. 

Kelbg, G.: Ann. Physik, 12 (1963) 219; 13 (1964) 354; 14 (1964) 394 
Ebeling, W., Richert, W., Kraeft, W. D., Stolzmann, W.: phys.stat.sol.(b) 104 (1981) 
193 

Haronska, P., Kremp, D., Schlanges, M: Wiss. Z. Univ. Rostock, MN Reihe 36 (1987) 
98 

Ebeling, W., Richert, W.: Phys. Lett. A 108 (1985) 80; phys. stat sol. (b) 128 (1985) 
67 



